BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161114 16.txt","14 Nov 2016 22:31",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161114 17.txt","14 Nov 2016 23:01",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161114 18.txt","14 Nov 2016 23:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 19.txt","15 Nov 2016 00:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 20.txt","15 Nov 2016 00:31",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 21.txt","15 Nov 2016 01:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 22.txt","15 Nov 2016 01:32",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 23.txt","15 Nov 2016 01:59",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 24.txt","15 Nov 2016 02:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 25.txt","15 Nov 2016 02:59",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 26.txt","15 Nov 2016 03:29",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 27.txt","15 Nov 2016 04:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 28.txt","15 Nov 2016 04:33",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 29.txt","15 Nov 2016 05:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 30.txt","15 Nov 2016 05:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 31.txt","15 Nov 2016 05:59",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 32.txt","15 Nov 2016 06:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 33.txt","15 Nov 2016 06:59",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 34.txt","15 Nov 2016 07:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 35.txt","15 Nov 2016 08:07",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 36.txt","15 Nov 2016 08:29",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 37.txt","15 Nov 2016 08:59",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 38.txt","15 Nov 2016 09:31",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 39.txt","15 Nov 2016 10:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 40.txt","15 Nov 2016 10:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 41.txt","15 Nov 2016 11:01",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 42.txt","15 Nov 2016 11:31",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 43.txt","15 Nov 2016 12:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 44.txt","15 Nov 2016 12:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 45.txt","15 Nov 2016 13:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 46.txt","15 Nov 2016 13:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 47.txt","15 Nov 2016 14:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 48.txt","15 Nov 2016 14:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 49.txt","15 Nov 2016 15:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 50.txt","15 Nov 2016 15:30",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov <- ReadProfiles("1_BRAM3_Inclinometer_20161115 51.txt","15 Nov 2016 16:00",FALSE,
BramleyNov$label,BramleyNov$Depth,BramleyNov$A0,BramleyNov$A180,BramleyNov$B0,BramleyNov$B180)
BramleyNov$Depth <- (BramleyNov$Depth - 0.5)   #Excluding height of pully wheel
save(BramleyNov,file="BramleyNov.RData")
load("BramleyNov.RData")
# Process inclinometer data
# Initialise the functions
source("displacement.R")
source("PlotProfiles.R")
source("surveydiff.R")
# Close existing plot windows if any are open
px <- length(dev.list())
if (px>0) {for (ix in 1:px) {dev.off()}}
# Constants
const <- 25000		# Multiplier used by DataMate to scale sine of angles
measure <-500 		# Distance between measurements (mm)
TY <-  c("15 Nov 2016 06:30","14 Nov 2016 14:32","14 Nov 2016 15:21","14 Nov 2016 15:43","14 Nov 2016 16:11","14 Nov 2016 16:54","14 Nov 2016 17:29","14 Nov 2016 17:44",
"14 Nov 2016 18:31","14 Nov 2016 19:31","14 Nov 2016 20:01","14 Nov 2016 20:31","14 Nov 2016 20:59","14 Nov 2016 21:34",
"14 Nov 2016 22:00","14 Nov 2016 22:31","14 Nov 2016 23:01","14 Nov 2016 23:30","15 Nov 2016 00:00","15 Nov 2016 00:31","15 Nov 2016 01:00",
"15 Nov 2016 01:32","15 Nov 2016 01:59","15 Nov 2016 02:30","15 Nov 2016 02:59","15 Nov 2016 03:29","15 Nov 2016 04:00","15 Nov 2016 04:33",
"15 Nov 2016 05:00","15 Nov 2016 05:30","15 Nov 2016 05:59","15 Nov 2016 06:30","15 Nov 2016 06:59","15 Nov 2016 07:30","15 Nov 2016 08:07",
"15 Nov 2016 08:29","15 Nov 2016 08:59","15 Nov 2016 09:31","15 Nov 2016 10:00","15 Nov 2016 10:30","15 Nov 2016 11:01","15 Nov 2016 11:31",
"15 Nov 2016 12:00","15 Nov 2016 12:30","15 Nov 2016 13:00","15 Nov 2016 13:30","15 Nov 2016 14:00","15 Nov 2016 14:30","15 Nov 2016 15:00",
"15 Nov 2016 15:30","15 Nov 2016 16:00")
TYDay <- strptime(TY,"%d %b %Y %H:%M", tz="NZ")
levels <- c(1,11,21,31)
# Load the profile data and assemble into data frames for analysis
# Tyler's continuous 12 hour data
load("BramleyNov.RData.RData")
# Process the data - first select dataset
Profiles <- BramleyNov
tDay <- TYDay
# How many profiles are present?
nx <- length(Profiles$label)
nn <- nx-1
nl <- length(levels)
ns <- length(tDay)
if (ns!=nx) {print(paste("Number of profiles,",nx,", does not equal number of measurement times,",ns))}
# Loop through all the profiles and calculate displacements
for (ix in 1:nx) {
if (ix==1) {
offsets <- displacement(Profiles$A0[,ix],Profiles$A180[,ix],Profiles$B0[,ix],Profiles$B180[,ix])
} else {
temp <- displacement(Profiles$A0[,ix],Profiles$A180[,ix],Profiles$B0[,ix],Profiles$B180[,ix])
offsets$Adiff <- cbind(offsets$Adiff,temp$Adiff)
offsets$Bdiff <- cbind(offsets$Bdiff,temp$Bdiff)
offsets$Adev <- cbind(offsets$Adev,temp$Adev)
offsets$Bdev <- cbind(offsets$Bdev,temp$Bdev)
offsets$Ack <- cbind(offsets$Ack,temp$Ack)
offsets$Bck <- cbind(offsets$Bck,temp$Bck)
offsets$Asum <- cbind(offsets$Asum,temp$Asum)
offsets$Bsum <- cbind(offsets$Bsum,temp$Bsum)
}
}
# Now calculate the displacements between surveys if nx>1 and plot results if plotit=TRUE
plotit = FALSE
if (nx>1) {
for (ix in 2:nx) {
if (ix==2) {
trend <- surveydiff(Profiles,offsets,ix,plotit)
} else {
temp <- surveydiff(Profiles,offsets,ix,plotit)
trend$Adisp <- cbind(trend$Adisp,temp$Adisp)
trend$Bdisp <- cbind(trend$Bdisp,temp$Bdisp)
trend$Asum <- cbind(trend$Asum,temp$Asum)
trend$Bsum <- cbind(trend$Bsum,temp$Bsum)
}
}
}
# Now extract slices
SliceA <- matrix(data=NA,nrow=nl,ncol=nn)
SliceB <- matrix(data=NA,nrow=nl,ncol=nn)
SliceI <- matrix(data=NA,nrow=nl,ncol=1)
for (ix in 1:nl) {
SliceI[ix] <- which(Profiles$Depth[,2]==levels[ix])
}
for (ix in 1:nn) {
for (iy in 1:nl) {
SliceA[iy,ix] <- trend$Asum[SliceI[iy],ix]
SliceB[iy,ix] <- trend$Bsum[SliceI[iy],ix]
}
}
# Plot profiles
strt <- 1
for (ix in strt:nx) {
PlotProfiles(Profiles,offsets,trend,ix)
}
# Plot time series of displacements
# As profiles
dev.new()
op <- par(mfcol=c(1,2))
xla <- max(abs(trend$Asum),2)
cols <- 1
ltys <- 1
plot(trend$Asum[,1],Profiles$Depth[,2],xlim=c(-xla,xla),ylim=c(43,0),type="l",col=1,lty=1,
xlab="Deviation (mm)",ylab="Depth (m)",
main="Displacement - A axis",sub = paste("Since",Profiles$label[1]))
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
lines(trend$Asum[,ix],Profiles$Depth[,(ix+1)],col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
legend("left",Profiles$label[2:nx],col=cols,lty=ltys,cex=0.8,bg="white")
xlb <- max(abs(trend$Bsum),2)
cols <- 1
ltys <- 1
plot(trend$Bsum[,1],Profiles$Depth[,(ix+1)],xlim=c(-xlb,xlb),ylim=c(43,0),type="l",col=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50),lty=1,
xlab="Deviation (mm)",ylab="Depth (m)",
main="Displacement - B axis",sub = paste("Since",Profiles$label[1]))
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
lines(trend$Bsum[,ix],Profiles$Depth[,(ix+1)],col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
#legend("bottomright",Profiles$label[2:nx],col=cols,lty=ltys,cex=0.75,bg="white")
par(op)
# As a plan view #showing plan view with all data
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)+2
dev.new()
cols <- 1
ltys <- 1
pchs <- 1
plot(trend$Bsum[,1],trend$Asum[,1],type="o",pch=1,col=1,lty=1,xlim=c(-xlc,xlc),ylim=c(-xlc,xlc),
main=paste("Displacement since",Profiles$label[1]),xlab="B axis (mm)",ylab="A axis (mm)")
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
pchs <- cbind(pchs,ix)
lines(trend$Bsum[,ix],trend$Asum[,ix],type="o",pch=ix,col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(h=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,lty=ltys,cex=0.75,bg="white")
# As slices # displacemnt over time split into A and B
dev.new()
op <- par(mfcol=c(2,1))
cols <- 1
pchs <- 1
ltys <- 1
plot(tDay[2:ns],SliceA[1,],type="o",pch=pchs,col=cols,lty=ltys,ylim=c(-2,2),
main="A axis displacements",xlab="Time",ylab="Displacement (mm)")
for (ix in 2:nl) {
cols <- cbind(cols,ix)
pchs <- cbind(pchs,ix)
ltys <- cbind(pchs,ix)
lines(tDay[2:ns],SliceA[ix,],type="o",pch=ix,col=ix,lty=ix)
}
abline(h=0,lty=2,col="black")
legend("bottomleft",as.character(levels),col=cols,pch=pchs,cex=0.5,bg="white")
plot(tDay[2:ns],SliceB[1,],type="o",pch=pchs[1],col=cols[1],lty=ltys[1],ylim=c(-3,6),
main="B axis displacements",xlab="Time",ylab="Displacement (mm)")
for (ix in 2:nl) {
points(tDay[2:ns],SliceB[ix,],type="o",pch=pchs[ix],col=pchs[ix],lty=ltys[ix])
}
abline(h=0,lty=2,col="black")
legend("topleft",as.character(levels),col=cols,pch=pchs,lty=ltys,cex=0.5,bg="white")
# plot(moon.distance$Distance~moon.distance$Time, type="o",xlim=as.POSIXct(c("2013-07-01 18:50","2016-05-15 14:47")))
par(op)
# As a plan view for 1 m level #displacment moving around over time
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)
ilevel <- which(Profiles$Depth[,1]==levels[1])
for (ix in 2:length(levels)) {
ilevel <- cbind(ilevel,which(Profiles$Depth[,1]==levels[ix]))
}
dev.new()
cols <- 1
pchs <- 1
plot(trend$Bsum[1,1],trend$Asum[1,1],type="o",pch=1,col=1,xlim=c(-2,4),ylim=c(-2,4),
main=paste("Displacement since",Profiles$label[1]," at ",levels[1]," m level"),
xlab="B axis (mm)",ylab="A axis (mm)")
xBsum <- trend$Bsum[1,1]
yAsum <- trend$Asum[1,1]
for (ix in 2:nn) {
cols <- cbind(cols,ix)
if (ix<26) pchs <- cbind(pchs,ix) else pchs <- cbind(pchs-25,ix)
lines(trend$Bsum[1,ix],trend$Asum[1,ix],type="o",pch=ix,col=ix)
xBsum <- cbind(xBsum,trend$Bsum[1,ix])
yAsum <- cbind(yAsum,trend$Asum[1,ix])
}
lines(xBsum,yAsum,lty=1,col="black")
abline(v=0,lty=2,col="dark grey")
abline(h=0,lty=2,col="dark grey")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,cex=0.75,bg="white")
# Process inclinometer data
# Initialise the functions
source("displacement.R")
source("PlotProfiles.R")
source("surveydiff.R")
# Close existing plot windows if any are open
px <- length(dev.list())
if (px>0) {for (ix in 1:px) {dev.off()}}
# Constants
const <- 25000		# Multiplier used by DataMate to scale sine of angles
measure <-500 		# Distance between measurements (mm)
TY <-  c("15 Nov 2016 06:30","14 Nov 2016 14:32","14 Nov 2016 15:21","14 Nov 2016 15:43","14 Nov 2016 16:11","14 Nov 2016 16:54","14 Nov 2016 17:29","14 Nov 2016 17:44",
"14 Nov 2016 18:31","14 Nov 2016 19:31","14 Nov 2016 20:01","14 Nov 2016 20:31","14 Nov 2016 20:59","14 Nov 2016 21:34",
"14 Nov 2016 22:00","14 Nov 2016 22:31","14 Nov 2016 23:01","14 Nov 2016 23:30","15 Nov 2016 00:00","15 Nov 2016 00:31","15 Nov 2016 01:00",
"15 Nov 2016 01:32","15 Nov 2016 01:59","15 Nov 2016 02:30","15 Nov 2016 02:59","15 Nov 2016 03:29","15 Nov 2016 04:00","15 Nov 2016 04:33",
"15 Nov 2016 05:00","15 Nov 2016 05:30","15 Nov 2016 05:59","15 Nov 2016 06:30","15 Nov 2016 06:59","15 Nov 2016 07:30","15 Nov 2016 08:07",
"15 Nov 2016 08:29","15 Nov 2016 08:59","15 Nov 2016 09:31","15 Nov 2016 10:00","15 Nov 2016 10:30","15 Nov 2016 11:01","15 Nov 2016 11:31",
"15 Nov 2016 12:00","15 Nov 2016 12:30","15 Nov 2016 13:00","15 Nov 2016 13:30","15 Nov 2016 14:00","15 Nov 2016 14:30","15 Nov 2016 15:00",
"15 Nov 2016 15:30","15 Nov 2016 16:00")
TYDay <- strptime(TY,"%d %b %Y %H:%M", tz="NZ")
levels <- c(1,11,21,31)
# Load the profile data and assemble into data frames for analysis
# Tyler's continuous 12 hour data
load("BramleyNov.RData")
# Process the data - first select dataset
Profiles <- BramleyNov
tDay <- TYDay
# How many profiles are present?
nx <- length(Profiles$label)
nn <- nx-1
nl <- length(levels)
ns <- length(tDay)
if (ns!=nx) {print(paste("Number of profiles,",nx,", does not equal number of measurement times,",ns))}
# Loop through all the profiles and calculate displacements
for (ix in 1:nx) {
if (ix==1) {
offsets <- displacement(Profiles$A0[,ix],Profiles$A180[,ix],Profiles$B0[,ix],Profiles$B180[,ix])
} else {
temp <- displacement(Profiles$A0[,ix],Profiles$A180[,ix],Profiles$B0[,ix],Profiles$B180[,ix])
offsets$Adiff <- cbind(offsets$Adiff,temp$Adiff)
offsets$Bdiff <- cbind(offsets$Bdiff,temp$Bdiff)
offsets$Adev <- cbind(offsets$Adev,temp$Adev)
offsets$Bdev <- cbind(offsets$Bdev,temp$Bdev)
offsets$Ack <- cbind(offsets$Ack,temp$Ack)
offsets$Bck <- cbind(offsets$Bck,temp$Bck)
offsets$Asum <- cbind(offsets$Asum,temp$Asum)
offsets$Bsum <- cbind(offsets$Bsum,temp$Bsum)
}
}
# Now calculate the displacements between surveys if nx>1 and plot results if plotit=TRUE
plotit = FALSE
if (nx>1) {
for (ix in 2:nx) {
if (ix==2) {
trend <- surveydiff(Profiles,offsets,ix,plotit)
} else {
temp <- surveydiff(Profiles,offsets,ix,plotit)
trend$Adisp <- cbind(trend$Adisp,temp$Adisp)
trend$Bdisp <- cbind(trend$Bdisp,temp$Bdisp)
trend$Asum <- cbind(trend$Asum,temp$Asum)
trend$Bsum <- cbind(trend$Bsum,temp$Bsum)
}
}
}
# Now extract slices
SliceA <- matrix(data=NA,nrow=nl,ncol=nn)
SliceB <- matrix(data=NA,nrow=nl,ncol=nn)
SliceI <- matrix(data=NA,nrow=nl,ncol=1)
for (ix in 1:nl) {
SliceI[ix] <- which(Profiles$Depth[,2]==levels[ix])
}
for (ix in 1:nn) {
for (iy in 1:nl) {
SliceA[iy,ix] <- trend$Asum[SliceI[iy],ix]
SliceB[iy,ix] <- trend$Bsum[SliceI[iy],ix]
}
}
# Plot profiles
strt <- 1
for (ix in strt:nx) {
PlotProfiles(Profiles,offsets,trend,ix)
}
# Plot time series of displacements
# As profiles
dev.new()
op <- par(mfcol=c(1,2))
xla <- max(abs(trend$Asum),2)
cols <- 1
ltys <- 1
plot(trend$Asum[,1],Profiles$Depth[,2],xlim=c(-xla,xla),ylim=c(43,0),type="l",col=1,lty=1,
xlab="Deviation (mm)",ylab="Depth (m)",
main="Displacement - A axis",sub = paste("Since",Profiles$label[1]))
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
lines(trend$Asum[,ix],Profiles$Depth[,(ix+1)],col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
legend("left",Profiles$label[2:nx],col=cols,lty=ltys,cex=0.8,bg="white")
xlb <- max(abs(trend$Bsum),2)
cols <- 1
ltys <- 1
plot(trend$Bsum[,1],Profiles$Depth[,(ix+1)],xlim=c(-xlb,xlb),ylim=c(43,0),type="l",col=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50),lty=1,
xlab="Deviation (mm)",ylab="Depth (m)",
main="Displacement - B axis",sub = paste("Since",Profiles$label[1]))
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
lines(trend$Bsum[,ix],Profiles$Depth[,(ix+1)],col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
#legend("bottomright",Profiles$label[2:nx],col=cols,lty=ltys,cex=0.75,bg="white")
par(op)
# As a plan view #showing plan view with all data
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)+2
dev.new()
cols <- 1
ltys <- 1
pchs <- 1
plot(trend$Bsum[,1],trend$Asum[,1],type="o",pch=1,col=1,lty=1,xlim=c(-xlc,xlc),ylim=c(-xlc,xlc),
main=paste("Displacement since",Profiles$label[1]),xlab="B axis (mm)",ylab="A axis (mm)")
for (ix in 2:nn) {
cols <- cbind(cols,ix)
ltys <- cbind(ltys,ix)
pchs <- cbind(pchs,ix)
lines(trend$Bsum[,ix],trend$Asum[,ix],type="o",pch=ix,col=ix,lty=ix)
}
abline(v=0,lty=2,col="black")
abline(h=0,lty=2,col="black")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,lty=ltys,cex=0.75,bg="white")
# As slices # displacemnt over time split into A and B
dev.new()
op <- par(mfcol=c(2,1))
cols <- 1
pchs <- 1
ltys <- 1
plot(tDay[2:ns],SliceA[1,],type="o",pch=pchs,col=cols,lty=ltys,ylim=c(-2,2),
main="A axis displacements",xlab="Time",ylab="Displacement (mm)")
for (ix in 2:nl) {
cols <- cbind(cols,ix)
pchs <- cbind(pchs,ix)
ltys <- cbind(pchs,ix)
lines(tDay[2:ns],SliceA[ix,],type="o",pch=ix,col=ix,lty=ix)
}
abline(h=0,lty=2,col="black")
legend("bottomleft",as.character(levels),col=cols,pch=pchs,cex=0.5,bg="white")
plot(tDay[2:ns],SliceB[1,],type="o",pch=pchs[1],col=cols[1],lty=ltys[1],ylim=c(-3,6),
main="B axis displacements",xlab="Time",ylab="Displacement (mm)")
for (ix in 2:nl) {
points(tDay[2:ns],SliceB[ix,],type="o",pch=pchs[ix],col=pchs[ix],lty=ltys[ix])
}
abline(h=0,lty=2,col="black")
legend("topleft",as.character(levels),col=cols,pch=pchs,lty=ltys,cex=0.5,bg="white")
# plot(moon.distance$Distance~moon.distance$Time, type="o",xlim=as.POSIXct(c("2013-07-01 18:50","2016-05-15 14:47")))
par(op)
# As a plan view for 1 m level #displacment moving around over time
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)
ilevel <- which(Profiles$Depth[,1]==levels[1])
for (ix in 2:length(levels)) {
ilevel <- cbind(ilevel,which(Profiles$Depth[,1]==levels[ix]))
}
dev.new()
cols <- 1
pchs <- 1
plot(trend$Bsum[1,1],trend$Asum[1,1],type="o",pch=1,col=1,xlim=c(-2,4),ylim=c(-2,4),
main=paste("Displacement since",Profiles$label[1]," at ",levels[1]," m level"),
xlab="B axis (mm)",ylab="A axis (mm)")
xBsum <- trend$Bsum[1,1]
yAsum <- trend$Asum[1,1]
for (ix in 2:nn) {
cols <- cbind(cols,ix)
if (ix<26) pchs <- cbind(pchs,ix) else pchs <- cbind(pchs-25,ix)
lines(trend$Bsum[1,ix],trend$Asum[1,ix],type="o",pch=ix,col=ix)
xBsum <- cbind(xBsum,trend$Bsum[1,ix])
yAsum <- cbind(yAsum,trend$Asum[1,ix])
}
lines(xBsum,yAsum,lty=1,col="black")
abline(v=0,lty=2,col="dark grey")
abline(h=0,lty=2,col="dark grey")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,cex=0.75,bg="white")
# As a plan view for 1 m level #displacment moving around over time
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)
ilevel <- which(Profiles$Depth[,1]==levels[1])
for (ix in 2:length(levels)) {
ilevel <- cbind(ilevel,which(Profiles$Depth[,1]==levels[ix]))
}
dev.new()
cols <- 1
pchs <- 1
plot(trend$Bsum[1,1],trend$Asum[1,1],type="o",pch=1,col=1,xlim=c(-4,2),ylim=c(-2,4),
main=paste("Displacement since",Profiles$label[1]," at ",levels[1]," m level"),
xlab="B axis (mm)",ylab="A axis (mm)")
xBsum <- trend$Bsum[1,1]
yAsum <- trend$Asum[1,1]
for (ix in 2:nn) {
cols <- cbind(cols,ix)
if (ix<26) pchs <- cbind(pchs,ix) else pchs <- cbind(pchs-25,ix)
lines(trend$Bsum[1,ix],trend$Asum[1,ix],type="o",pch=ix,col=ix)
xBsum <- cbind(xBsum,trend$Bsum[1,ix])
yAsum <- cbind(yAsum,trend$Asum[1,ix])
}
lines(xBsum,yAsum,lty=1,col="black")
abline(v=0,lty=2,col="dark grey")
abline(h=0,lty=2,col="dark grey")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,cex=0.75,bg="white")
# As a plan view for 1 m level #displacment moving around over time
xlc <- max(abs(trend$Asum),abs(trend$Bsum),2)
ilevel <- which(Profiles$Depth[,1]==levels[1])
for (ix in 2:length(levels)) {
ilevel <- cbind(ilevel,which(Profiles$Depth[,1]==levels[ix]))
}
dev.new()
cols <- 1
pchs <- 1
plot(trend$Bsum[1,1],trend$Asum[1,1],type="o",pch=1,col=1,xlim=c(-5,2),ylim=c(-2,4),
main=paste("Displacement since",Profiles$label[1]," at ",levels[1]," m level"),
xlab="B axis (mm)",ylab="A axis (mm)")
xBsum <- trend$Bsum[1,1]
yAsum <- trend$Asum[1,1]
for (ix in 2:nn) {
cols <- cbind(cols,ix)
if (ix<26) pchs <- cbind(pchs,ix) else pchs <- cbind(pchs-25,ix)
lines(trend$Bsum[1,ix],trend$Asum[1,ix],type="o",pch=ix,col=ix)
xBsum <- cbind(xBsum,trend$Bsum[1,ix])
yAsum <- cbind(yAsum,trend$Asum[1,ix])
}
lines(xBsum,yAsum,lty=1,col="black")
abline(v=0,lty=2,col="dark grey")
abline(h=0,lty=2,col="dark grey")
abline(v=c(-1.4,1.4),lty=4,col="grey")
abline(h=c(-1.4,1.4),lty=4,col="grey")
legend("bottomright",Profiles$label[2:nx],col=cols,pch=pchs,cex=0.75,bg="white")
###moon distance to be used in other scripts
library(zoo)
moon.distance <- read.csv("moon distance.csv", sep=";")
colnames(moon.distance) <- c("Time","Distance") # changing coloums names
moon.distance$Time <- as.POSIXct(moon.distance$Time, format="%Y-%b-%d %H:%M", tz="UTC")
moon.distance$Time <- format(moon.distance$Time, tz="NZ")
moon.distance$Time <- as.POSIXct(moon.distance$Time, tz="NZ")
moon.distance <- na.omit(moon.distance)
moon.distance.zoo <- zoo(moon.distance[,2], moon.distance[,1])
plot(moon.distance.zoo)
